import load_foil as lf
import numpy as np
import matplotlib.pyplot as mp

mp.clf()

#-------chargement de l aile sous forme d un vecteur de points en 2D----------#
(ex,ey,ix,iy)= lf.load_foil("BACNLF.DAT")
mp.plot(ex,ey)
mp.plot(ix,iy)
mp.title("Wing curve without Cubic Spline")
mp.xlabel("Abscissa x")
mp.ylabel("Ordonate y")
mp.savefig("Sans-spline.png")

#-------recuperation de la taille des vecteurs representant respectivement la partie sup et inf de l aile----#
ne = ex.shape[0]
ni = ix.shape[0]

#------choix aleatoire des valeurs de la derivee de la fonction y aux points x1 et xn pour
#------permettre la resolution du systeme lineaire de y"---------------------------------#
yp1=0.0
ypn=0.0




#-------Fonction retournant la table des derivees secondes
#-------de y en resolvant le systeme triangulaire--------#
#O(n)
def spline(ex,ey,n,yp1,ypn):
    u=np.zeros([1,n-1])
    u=u[0]
    y2 = np.zeros([1,n])
    y2 = y2[0]
    if(yp1>0.99e30):
        y2[0]=0.0
        u[0]=0.0
     
     
    else :
        y2[0] = -0.5
        u[0]=(3.0/(ex[2]-ex[1]))*((ey[2]-ey[1])/(ex[2]-ex[1])-yp1)
        
        
    for i in np.arange(1,n-1,1):
        sig=(ex[i]-ex[i-1])/(ex[i+1]-ex[i-1])
        p=sig*y2[i-1]+2.0
        y2[i]=(sig-1.0)/p
        u[i]=(ey[i+1]-ey[i])/(ex[i+1]-ex[i])-(ey[i]-ey[i-1])/(ex[i]-ex[i-1])
        u[i]=(6.0*u[i]/(ex[i+1]-ex[i-1])-sig*u[i-1])/p
        
        
    if (ypn>0.99e30):
        qn=0.0
        un=0.0
    else :
        qn = 0.5
        un =(3.0/(ex[n-1]-ex[n-2]))*(ypn-(ey[n-1]-ey[n-2])/(ex[n-1]-ex[n-2]))
       
        
    y2[n-1]=(un-qn*u[n-2])/(qn*y2[n-2]+1.0)
    for k in np.arange(n-2,-1,-1):
        y2[k]=y2[k]*y2[k+1]+u[k]

    return y2




#-------Fonction calculant la nouvelle valeur de la fonction y apres interpolation
#----------------------------a un point x donne----------------------------------#
#O(log2(n))
def splint(xa, ya , n, x , y2):
   
    klo = 0.0
    khi = n-1
    k = 0

    while(khi - klo > 1):
        k = (khi + klo) / 2
        if(xa[k] > x):
            khi=k
        else:
            klo=k
    h = xa[khi]-xa[klo]
    if(h == 0.0):
        print "Bad xa input to routine splint"
    a = (xa[khi] - x) / h
    b = (x - xa[klo]) / h
    
    y = a*ya[klo] + b*ya[khi] + ((a**3 - a)*y2[klo] + (b**3 - b)*y2[khi])*(h**2)/6.0

    return y


#-----Fonction permettant de constituer les vecteurs des differents valeurs de la fonction y
#-----en un nombre determine de points figurant dans  un interval donne--------------------#
#O(N)
def cubic_spline(xa , ya , n , y2 , N ,i):
    X = np.asarray(np.zeros([N,1]))
    Y = np.asarray(np.zeros([N,1]))
    pas = (xa[i+1]-xa[i])/(N-1)
    X[0] = xa[i]
    Y[0] = ya[i]
    k = 1
    for z in np.arange(xa[i]+pas,xa[i+1],pas):
        X[k] = z
        Y[k] = splint(xa, ya, n, z , y2)
        k = k+1
   
    
    X[N-1] = xa[i+1]
    Y[N-1] = ya[i+1]
    
    return (X,Y)
    
    
#------Fonction permettant de tracer la courbe d interpolation interval par interval----#
def draw_spline_cubic(xa , ya , n , N , y2):
    
    for i in range(n-1):
        (X,Y) = cubic_spline(xa , ya , n , y2 , N ,i)
        mp.plot(X , Y, 'r' )
        mp.plot(np.array([X[0]]) , np.array([Y[0]]) , 'b' ,markeredgewidth=1,markersize=5,marker='o')
        
mp.clf()
draw_spline_cubic(ex,ey,ne,10, spline(ex,ey,ne,yp1,ypn))
draw_spline_cubic(ix,iy,ni,10, spline(ix,iy,ni,yp1,ypn))
mp.title("Wing curve using Cubic Spline with y'(x1)=y'(xn)=0.0")
mp.xlabel("Abscissa x")
mp.ylabel("Ordonate y")

mp.savefig("Avec-spline1.png")



mp.clf()
draw_spline_cubic(ex,ey,ne,10, spline(ex,ey,ne,7,7.4))
draw_spline_cubic(ix,iy,ni,10, spline(ix,iy,ni,7,7.4))
mp.title("Wing curve using Cubic Spline with y'(x1)=7 and y'(xn)=7.4")
mp.xlabel("Abscissa x")
mp.ylabel("Ordonate y")

mp.savefig("Avec-spline2.png")
















